      subroutine refprofs
c-----------------------------------------------------------------------
c  specify vertical profile of reference variables
c  this subrout interpolates a profile. To speed
c  up the process the data is interpolated to a regularly spaced
c  grid. to avoid errors, all the points in the given profile should
c  belong to the regular grid.
c
c ioreft: 1 - interpolate thdat as pts;
c         2 - interpolate thdat as tems;
c         3 - interpolate thdat as N**2 and calculates pts profile.
c
c iorefq  1 - interpolate qvdat as relative humidity (0..1);
c         2 - interpolate qvdat as specific humidity (kg/kg);
c
c ptdat is a vertical profile of potential temperature calculated
c in the case ioreft.ne.2 to compute the pressure profile.
c
c note: the interpolation is done using z as the vertical coordinate.
c in fact, reference profiles should be a function of p.
c-----------------------------------------------------------------------
c p.m. 2005/10/11
c-----------------------------------------------------------------------

      use alloc
      use refstate
      implicit none
      integer i,ip,id,iz
      double precision bvf2m,ptmed,qsat,tid,tmed
      character*80 fndebug
      integer istepp
      save istepp
      data istepp/0/

c initialization of the profile: zdat is (in principle) irregularly
c spaced; zpro is regularly spaced. this should be done only once.

      if(verbose.ge.2) write(*,*) 'Refprofs ',npro,dzpro,nprof,ndat
     :  ,pts0
      if(istepp.eq.0) then
        do i=0,npro
          zpro(i)=i*dzpro
        enddo

c compute weighting factors for profiles

        fndebug='weiprof.grd'
        call weiprof(fndebug)
        istepp=1
      endif

      do ip=1,nprof

c NEW: the surface pressure is read from the ref file

      psref=pressdat(0,ip)
      press0(ip)=pressdat(0,ip)
      if(ioreft.eq.1) then
        call inipro(thdat(0,ip),zthdat(0,ip),ndat,fpro(0),zpro,npro)
        t0ini=thdat(0,ip)*(psref/p00)**akapa
        if(verbose.ge.2) write(*,*) 't0ini=',t0ini
        do id=0,ndat
          ptdat(id,ip)=thdat(id,ip)
        enddo
      elseif(ioreft.eq.2) then
        call inipro(thdat(0,ip),zthdat(0,ip),ndat,fpro(0),zpro,npro)
        t0ini=thdat(0,ip)
        if(verbose.ge.3) then
          do i=0,npro
            write(*,*) 'th:',i,fpro(i)
          enddo
        endif
        do id=0,ndat
          ptdat(id,ip)=thdat(id,ip)*(pressdat(id,ip)/1.e5)**(-akapa)
        enddo
      else
        call inipro(thdat(0,ip),zthdat(0,ip),ndat,fnpro(0),zpro,npro)
        t0ini=pts0*(psref/p00)**akapa
        if(verbose.ge.2) write(*,*) 't0ini=',t0ini,fnpro(1)
c
        fpro(0)=pts0
        do iz=1,npro
          bvf2m=0.5*(fnpro(iz-1)+fnpro(iz))
          fpro(iz)=fpro(iz-1)*dexp(bvf2m*dzpro/g)
c         if(verbose.ge.2) write(*,*) 'th:',iz,bvf2m,fpro(iz)
        enddo
        fpro(npro+1)=fpro(npro)
c
        ptdat(0,ip)=t0ini
        do id=1,ndat
          bvf2m=0.5*(thdat(id-1,ip)+thdat(id,ip))
          ptdat(id,ip)=ptdat(id-1,ip)
     :      *dexp(bvf2m*(zthdat(id,ip)-zthdat(id-1,ip))/g)
        enddo
      endif

      do i=0,npro
        dthdz(i,ip)=(fpro(i+1)-fpro(i))/dzpro
      enddo
c
      do i=0,npro
        thpro0(i,ip)=fpro(i)-dthdz(i,ip)*zpro(i)
      enddo

        if(verbose.ge.3) then
          do i=0,npro
            write(*,*) 'th:',i,thpro0(i,ip),dthdz(i,ip),zpro(i)
          enddo
        endif
c
      call uprof(iorefu,usdat(0,ip),zusdat(0,ip),ndat,fpro(0),zpro,npro)
c
      do i=0,npro
        dusdz(i,ip)=(fpro(i+1)-fpro(i))/dzpro
      enddo
c
      do i=0,npro
        uspro0(i,ip)=fpro(i)-dusdz(i,ip)*zpro(i)
      enddo
c
      call uprof(iorefv,vsdat(0,ip),zvsdat(0,ip),ndat,fpro(0),zpro,npro)
c
      do i=0,npro
        dvsdz(i,ip)=(fpro(i+1)-fpro(i))/dzpro
      enddo
c
      do i=0,npro
        vspro0(i,ip)=fpro(i)-dvsdz(i,ip)*zpro(i)
      enddo

c compute the pressure profile (integrate hydrostatic relation):

      psdat(0,ip)=psref
      if(verbose.ge.1) write(*,*) 'psref=',psref
      if(ioreft.eq.2) then

        if(iorefq.eq.1) then
          qvsdat(0,ip)=qvsdat(0,ip)*qsat(thdat(0,ip),psdat(0,ip))
        endif
        do id=1,ndat
          tmed=0.5*(thdat(id,ip)+thdat(id-1,ip))
          psdat(id,ip)=psdat(id-1,ip)*exp(-g/(r*tmed)
     :      *(zthdat(id,ip)-zthdat(id-1,ip)))

          if(iorefq.eq.1) then
            qvsdat(id,ip)=qvsdat(id,ip)
     :        *qsat(thdat(id,ip),psdat(id,ip))
          endif
          if(verbose.ge.2) write(*,*)
     :      id,' psdat=',psdat(id,ip),tmed,qvsdat(id,ip)
        enddo
      else
        tid=ptdat(0,ip)*(psdat(0,ip)/1.e5)**akapa
        if(iorefq.eq.1) then
          qvsdat(0,ip)=qvsdat(0,ip)*qsat(tid,psdat(0,ip))
        endif
        do id=1,ndat
          ptmed=0.5*(ptdat(id,ip)+ptdat(id-1,ip))
          psdat(id,ip)=(psdat(id-1,ip)**akapa-g/cp*p00**akapa
     :       *(zthdat(id,ip)-zthdat(id-1,ip))/ptmed)**(1./akapa)
          tid=ptdat(id,ip)*(psdat(id,ip)/1.e5)**akapa
          if(iorefq.eq.1) then
            qvsdat(id,ip)=qvsdat(id,ip)*qsat(tid,psdat(id,ip))
          endif
          if(verbose.ge.2) then
            write(*,*) id,' psdat=',psdat(id,ip),ptmed,qvsdat(id,ip)
     :        ,tid
          endif
        enddo
      endif

! compute temperature profile (needed for bc flux correction nhad0613ff)

      do id=0,ndat
        tedat(id,ip)=ptdat(id,ip)*(psdat(id,ip)/1.e5)**(akapa)
      enddo

      call inipro(tedat(0,ip),zthdat(0,ip),ndat,tepro(0,ip),zpro,npro)

      presspro(0,ip)=psref
      do i=1,npro
        presspro(i,ip)=presspro(i-1,ip)*exp(-g*(zpro(i)-zpro(i-1))
     :    /(r*0.5*(tepro(i,ip)+tepro(i-1,ip))))
      enddo

      do i=0,npro
        rhopro(i,ip)=presspro(i,ip)/(r*tepro(i,ip))
      enddo


      call inipro(qvsdat(0,ip),zqvsdat(0,ip),ndat,fpro(0),zpro,npro)
c
      do i=0,npro
        dqvsdz(i,ip)=(fpro(i+1)-fpro(i))/dzpro
      enddo
c
      do i=0,npro
        qvspro0(i,ip)=fpro(i)-dqvsdz(i,ip)*zpro(i)
      enddo
c
      write(nchn1,1000) ioreft,pts0,iorefu,iorefv,psref,ip
      do i=0,ndat
        write(nchn1,1001) zthdat(i,ip),thdat(i,ip)
     :    ,zusdat(i,ip),usdat(i,ip)
     :    ,zvsdat(i,ip),vsdat(i,ip),ptdat(i,ip)
     :    ,zqvsdat(i,ip),qvsdat(i,ip)
      enddo

      if(verbose.ge.2) then
        write(*,*) 'Interpolated profile ',ip
        write(*,'(9a15)') 'zpro','th','us','vs','qvs','dthdz'
     :    ,'temp','press','rho'
        do i=0,npro
          write(*,'(i5,9e15.7)')
     :      i,zpro(i),thpro0(i,ip),uspro0(i,ip),vspro0(i,ip)
     :        ,qvspro0(i,ip),dthdz(i,ip)
     :        ,tepro(i,ip),presspro(i,ip),rhopro(i,ip)
        enddo
      endif

      enddo
c
      return
1000  format(1x,'reference profiles: ioreft=',i2,' pts0=',f8.2
     :   ,' iorefu=',i2,' iorefv=',i2,' psref=',f7.0,' prof=',i3//
     :   ,t2,'zthdat',t10,'thdat',t25,'zusdat',t35,'usdat'
     :   ,t50,'zvsdat',t60,'vsdat',t75,'ptdat'//)
1001  format(1x,t2,f8.0,t10,e12.5,t25,f8.0,t35,f10.3,t50,f8.0,t60,f10.3
     :   ,t75,f10.3,t90,f10.3,t112,f10.5)
      end